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ABSTRACT 

A parallel algorithm for the efficient solution of a time dependent reaction convec- 
tion diffusion equation with small parameter on the diffusion term wdl be present^. 
The method is based on a domain decompos.t.on that is dictated by singular ? er 
turbation analysis. The analysis is used to determine regions where certain reduced 
= may be solved in place of the full equation. Parallelism s evident at two 
levels Domain decomposition provides parallelism at the highest level and within 
S domain there is ample opportunity to exploit parallelism. Run-t,me results 

demonstrate the viability of the method. 
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1. INTRODUCTION, In this paper, a new approach to solving partial differ- 
ential equations which model fluid flow is discussed and demonstrated. The algorithm 
is appropriate for modeling laminar transonic flow, such as through a duct of variable 
width. The method is an asymptotics-induced numerical method suitable for parallel 
processors which represent the state of the art in scientific computers. The contents 
of this paper concentrate on a description of the method and computational results. 
The complete theoretical basis for the algorithm has been developed in [32] and will 
appear separately. 

Competition between convection, diffusion, and reaction is crucial to the under- 
standing of fluid flow. When modeling transonic flow, except in regions of rapid 
variation such as in shocks and boundary layers, convection and reaction dominate 
over diffusion. A novel aspect of this method is the use of asymptotic analysis to 
exploit these physical properties, providing the theoretical basis for a domain decom- 
position. The analysis identifies the following two types of subdomains: regions where 
the solution is smooth, where a reduced equation may be solved; and regions of rapid 
variations, such as in a neighborhood of a shock, where the full equation must be 
solved. Domain decomposition provides large-grain parallelism. The domain decom- 
position is independent of the choice of numerical schemes for the subdomains; thus, 
schemes may be chosen which are a source of smaller-grain parallelism. Even though 
large grain parallelism is not exploited in the implementation, significant speedups 
are demonstrated. In addition to dictating the domain decomposition, asymptotics 
also provides a means of approximating solutions to the problem. In this way, a set of 
simplified problems is obtained that is better conditioned for numerical computations; 
hence, they may be solved by conventional techniques. The use of asymptotic analysis 
to precondition the computations is a new aspect of this method. 

The techniques presented herein are applicable to Computational Fluid Dynamics 
(CFD) in the transonic and supersonic regimes, in physical settings such as laminar 
flow through a nozzle (duct) and laminar flow around airfoils and other bodies. The 
gasdynamic equations, including viscous effects, are used as a model in these settings. 
Except for very simple geometries and boundary conditions there is no analytic solu- 
tion to these gasdynamic equations, and a numerical solution is difficult to obtain. For 
these reasons new algorithms are usually developed and tested on a more tractable 
canonical equation. The convection-diffusion-reaction equation 

du , . du d 2 u , v 

( 1 ) m +Mx ' t ’ U) te- ( d*- rix)u = 0 ' 

is such a canonical equation and will be the focus of this paper. The flows considered 
in this paper are not reacting fluids. Here, the reaction term arises from the effects of a 
variable cross sectional area in a duct. When the equation is nondimensionalized [24], 
the diffusion coefficient e is inversely proportional to the Reynolds number. Based 
on free-stream conditions in transonic flow, the Reynolds number for this problem 
is large. Asymptotic analysis exploits the smallness of the positive parameter e and 
involves study of the solution as c tends to zero (e i 0). This equation contains many 
of the properties that make the gasdynamic equations difficult to solve; namely, it is 
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capable of modeling rapid variations such as shocks and boundary layers. The method 
is capable of obtaining solutions to (1) when the shock is not stationary, which extends 
Howes studies [12,13] into the time-dependent regime. 

Asymptotic analysis gives qualitative and quantitative information as c j 0. The 
numerical method presented here exploits the analysis to determine an accurate so- 
lution for small positive e. The method is in the spirit of matched asymptotic ex- 
pansions [19,16], but it is not a numerical implementation of matched asymptotics. 
The asymptotic analysis involves the derivation of analytic upper and lower bounds 
on the solution, and is performed in the style of Howes [9,10,11]. Initially, bounds are 
discussed which are valid only in certain subregions. Then the bounds are combined 
to form a global a priori error bound. 

Another novel feature of the method is the availability of extensive error informa- 
tion in the form of both a priori error bounds and reliable a posteriori error estimates. 
Reliable a posteriori error estimates are obtained using the error analysis which ac- 
companies the numerical schemes used to solve the sub-problems. In addition, a priori 
error bounds are provided through the use of asymptotic analysis. The error analysis 
is based on the physical mechanisms associated with the problem; hence, it is based 
on accurate information (see [23]), not on the truncation of a Taylor series of a poorly 
behaved function. 

The method is an iterative technique. A linearized version of the original problem 
is solved in each step of the iteration. Theorems establishing the convergence of the 
method are presented, the proofs will appear in a subsequent paper. Computational 
experiments show that in just a few steps of the iteration, the solution to the nonlinear 
equation may be obtained. The iterative algorithm as well as the theorems associated 
with it are novel. 

In the next section, some of the ideas behind multiple scales asymptotic analysis 
are discussed. In addition, an introduction into how the asymptotic analysis and 
the numerical analysis are blended to form a computational method is presented. 
In Section 3 the problem is presented. Asymptotic analysis specific to this problem 
is discussed with the theorems supporting the method in Section 4. The iteration 
and method for detection of the subdomain boundary is discussed in Section 5. The 
numerical schemes used in the method are presented in Section 6. The method is 
stated in algorithmic form in Section 7. In Section 8 computational results on an 
Alliant FX/8 are presented. 


2. MULTIPLE SCALES. Many problems of scientific interest have multiple 
scales. These problems are characterized by the presence of distinguishable physical 
mechanisms, each associated with a temporal or spatial gauge or scale. When mod- 
eling a shock in a duct, for example, the width of the duct provides one scale, and 
the thickness of the shock layer provides another. The resolution of these scales is 
frequently required to determine the physics of interest. Asymptotic analysis provides 
analytic tools to identify and utilize the multiple scales. The relative importance of 
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any two physical processes in a given domain may be measured by the ratio of the 
corresponding scales; thus, the various scales may be ranked by a set of dimensionless 
parameters, the ratios of scales. When the ratio of two scales is a large or a small 
number, then it often happens that one of the competing mechanisms is dominant 
in most of the domain. For example, in laminar duct flow with large Reynolds num- 
ber the effects of viscosity may be ignored except in a neighborhood of the shock 
and boundary layers. The scales of the various competing processes (and, therefore, 
the relative magnitudes of the dimensionless parameters) usually change as the phe- 
nomenon evolves. Consider the behavior of the solution of the nonlinear parabolic 

equation, 

(2) P(u] := u t + uu x - eu xx - ru = 0, 

where e is a small positive parameter. This equation may be used as a model for shocks 
and boundary layers. For example, if r(x) = —a'(x)/a(x), where a(x) is the width of 
the duct of Figure 1, then this equation is associated with the flow through the duct 
[13]. There are (at least) two sets of scales appropriate when modeling shocks the 


A 
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FlG. 1. Variable width duct. (From [32].) 

scales associated with the original variables {x,t), and the scales appropriate in a 
small neighborhood of the shock (these are discussed in Section 4). 

The most easily tractable multiple-scale problems are those in which there are 
only a small number of widely separated groups of scales and the motion on the 
fastest scales has little influence on the smooth part of the solution. An identifying 
feature of this class is the presence of local regions in which the solution undergoes 
rapid variation. Such regions are called boundary or internal layers, when located in 
the neighborhood of a boundary or in the interior of the domain, respectively. These 
are the problems that are most natural for multitasking because it is easy to break up 
the domain according to the regions of different local behavior. The method presented 
here is appropriate for this class of multiple-scale problems. 

The decomposition into domains is accomplished using a symbiosis of numerics 
and asymptotics. The asymptotic analysis identifies the regions where diffusion is 
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negligible. In these regions, it is sufficient to solve a reduced equation. Solving this 
reduced equation can significantly reduce the work in the numerical method, and/or 
increase the potential for parallelism. For example, this allows the use of the method 
of characteristics to obtain a good approximation for the solution of (2). The numerics 
provides a means of solution in the subdomains, and also a feedback mechanism. The 
numerical scheme can expose regions of unexpected behavior, confirming or correcting 
the asymptotics-induced subdomain boundaries. This decomposition permits the use 
of locally refined meshes, allowing the concentration of computational effort in the 
regions where it is needed most. 

There is much literature on multiple-scales problems. Analytic methods for 
multiple-scales problems are discussed in the books [3,16,20,34]. The theory of multiple- 
scales analysis is discussed in [17,6,21]. The books [2,25] discuss both the techniques 
and the theory behind them. Finally, numerical techniques for multiple-scales prob- 
lems are discussed in the papers [22,26,33]. This list is meant only as an introduction 
to the literature, and not as a complete list. 


THE NONLINEAR PROBLEM. The method will be described and demon- 
strated by solving (2) on the domain 

( 3 ) D := {(x,f)|0 < x < 6,0 < t < T}, 
subject to 

( 4 ) u(x,0) = 7(x), 0 < x < b] 

( 5 ) u(0,f) = a(t), 0 < t < T; and 

( 6 ) u(b,t) = /3(t), 0<tCT. 

The portion of the boundary along which the data is specified is denoted by 
n := {(x,t)|0 < x < b, t = 0}f|{(a:,0|0<t < T, x = 0,6}. 

For the sake of simplicity, it is assumed that all boundaries are inflow boundaries, 

that is, a{t) > a 0 > 0 and 0(t) < 0 O < 0. The boundary data is assumed to be 
compatible; thus, 

( 7 ) Q (°)= 7(0 ), and 7(b) = /?(0). 

The coefficient of the forcing term r(x) is bounded with bounded derivatives. In 
addition, it is assumed that the solution to the reduced equation 

(8) P 0 [fi] := i it + &ii x — rh. = 0 
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has continuous derivatives at ( x,t ) = (0,0), and (x,t) — (6,0). This last restriction 
prevents the formation of corner layers, and may be expressed as 

(9) ^ + ^2x ~ n = °* f ° r ^ = ^°’ °^ ; 

(10) ^ + 7 £ _r7 = 0, for ( x » 4 ) = ( b ’°)' 

Under these conditions, the solution to (2) is uniquely defined [4]. 


4. ASYMPTOTIC ANALYSIS. Asymptotic analysis is employed to identify 
the dominant physics, creating an efficient and accurate numerical method. Here we 
sketch the analysis in the neighborhood of a shock. Readers wishing full detail and 
proofs of the results are encouraged to consult [32,11,10]. Shocks form in regions of 
merging characteristics (see Figure 2). Since the boundary conditions imposed on 



FIG. 2. Characteristics of steady-state solution u = — tanh[(.5 — x)/2e\, for e — .05. (From [32].) 

the problem are inflow conditions on both the x — 0 and x = b boundaries, the 
characteristics are traveling in the direction of increasing x from x = 0, and in the 
direction of decreasing x from x = b. These will merge (become asymptotically close) 
somewhere inside D , forming a shock. The merging of the characteristics stabilizes 
the shock, and keeps it from dispersing. 

Since the behavior of u as e | 0 is of interest, it is natural to first study the solution 
of the reduced equation (8). Weak solutions ft are sought for (8) with boundary data 

5 




(4-6). In order that & be uniquely defined, it is necessary to impose an entropy 
condition [18]. Suppose that it has a single shock. That is, suppose ft is the solution 
to ( 8 ) subject to (4-6) that is discontinuous only along a curve (x, t) = (r(i),f). For 
small c, this curve lies in the shock-layer region of the solution to the full problem. 
The size of this region tends to zero as e | 0 . Analytic methods for choosing T are 
discussed by Whitham [36], Kevorkian and Cole [16], and others. The path of the 
discontinuity is an analytic tool needed only for the theory. Since T is not needed for 
the computations, methods for choosing T will not be discussed here. 

The initial and boundary data are assumed to be smooth; thus, the shock does 
not exist at t = 0 . Rather, T is assumed to be undefined for t < t T , where t = t r is 
the time tz becomes discontinuous. It is natural to describe H in terms of the following 
functions: 

! Uo{x, t) for 0 < t < t r 

&i(x,t) for x < T - 1 (t) and t > t r 
fc 2 (x,f) forx>T _ 1 (t) and t > t r . 

The shock speed for u is 

(n) (s,(r(t),0 + s.(T(0,»))/*), 

so the entropy condition may be expressed as 

(12) Mr(t),<) > (Mr(t),t) + fc 2 (r(t),0)/2) > a 2 (r(f),t). 

Under these conditions, 

( 13 ) MO = *i(r(*),0 - fi 2 (r(0,t) > 0. 

The regions where fc is a good approximation to u are defined by presenting 
functions which bound the difference & — u. These bounds are small except in an 
asymptotically small neighborhood of the shock. The bounds, based on Howes [11], 
are reflected in the following theorem. 

THEOREM 1. (Howes) Let it be the solution to Po[&] = 0 and u the solution 
to P[u] = 0 on D, each satisfying the the boundary data ( 4 - 6 ). The solution to the 
reduced equation, is possibly discontinuous on the curve (x,f) = (r,t). Assume 
that the boundary data satisfy the following smoothness conditions: the data (4-6) 
satisfy the compatibility conditions (7), (9-10), and a, 0, 7 , with their first and second 
derivatives are all bounded. Then for e small enough 

|u tt | = 0(p exp[— / 2 (x, 0 / c ^ 2 ]) + 0( e) 

when the derivatives of i are continuous across T, and 

l u - &l = °{v exp[-/ 2 /c 1/2 ]) + 0(e 1/2 6 exp [-//e 1/2 ]) + O(e) 
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in the more general case when the derivatives of ft are not continuous across T. Here 
f{x,t) is a distance function between (x,f) and (I\t), and 6 is an upper bound on the 
difference of the normal derivative of it across T. 

It is now reasonable to utilize the theorem to make the definitions of the subdo- 
mains more precise. The internal layer is the following neighborhood of T: 

( 14 ) d il = {(x,OIOm) e D.lx-r- 1 ^)! < A(f)}. 

Here A(t) < Kr) (£) In 1 / 2 e is the width of the internal layer at time t ( K is a 
constant independent of e). The outer region is the complement of Du with respect 
to D, that is, 

(15) Dor = {(*, »)!(*,*) 6 A \x- > A(i)}. 

The upper bound on the size of the internal layer is based on the exp[-/ ! /V ,! ] terrn 
in the error bounds of Theorem 1. 

Theorem 1 motivates a preconditioning for the problem in Dor ■ The theorem 
states that (8) may be solved in place of (2). In addition, the theorem provides an 
error bound if diffusion (artificial or implicit in the numerical scheme) is incorporated 
into the solution process of either (8) or (2). Thus, the numerical method for D or 
may be chosen from the wide variety of methods designed for hyperbolic equations 

[1,35,8,28]. ,. 

The solution in the outer region is used to provide boundary data for the problem 

in the internal layer. (This is justified in Section 5.1.) Thus, it is possible to have 
boundary data for the internal-layer problem which is perturbed from the exact so- 
lution. The effects of this perturbation are that the height and location of the shock 
can vary by the same magnitude as the perturbation itself. This results in an error of 
magnitude 0(1) as e [ 0 in an asymptotically small neighborhood of the shock, and 
is reflected in the error bound of the following theorem [32]. 

THEOREM 2. Let u and v be solutions to (2) on the domain D IL with their 
boundary values satisfying the following smoothness conditions: the data arc bounded 
with bounded derivatives. Assume that the curve defining dD IL is smooth. Let 

i p(x,t) = u-v, for (x, t) 6 dDi L . 

Then for e and A small enough 

(16) \u - v\ = 0(ip) + 0((1 + f 2 /e 2 ) *), 

where f{x,t) measures the distance from ( x,t ) to (I\f). 

The theorem is the source for a second preconditioner. Namely, a scaled and 
translated coordinate system based on (16). Let / = jx-r(t)|. Setting x = (x- )/e, 

1 The curve is continuous with a continuous tangent. 
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the second term in the right hand side of (16) will be large when x = 0(1) as e | 0. An 
analogue of x will be used for the spatial coordinate in D IL , and will be described in 
Section 6.2. The use of this scaled and translated coordinate system creates a better 
conditioned numerical problem, thus it is a preconditioning. 

The local error bound of Theorem 1 is now used to form a global a priori error 
bound. The bound, as presented in Theorem 3 below, is sharp in D 0 r\ however, the 
bound reflects the possibility of shock displacement in the internal-layer region. 

THEOREM 3. Let u be the solution to (2) satisfying (4-6). Suppose v is obtained 
by first solving (8) in D 0 r subject to (4-6), then solving (2) on D IL with boundary 
data v on dD IL . Assume the compatibility conditions (7), (9-10) obtain, and that the 
data (4-6) are bounded with bounded derivates. If E = ||u - t;||i, then for e small 
enough 


E = O(e) 


in D or , and 


E = 0(e^ 4 In 1 / 2 c) 

in D il . 

The computational results were much stronger than the theorem suggests. Both 
the magnitude of A, and the magnitude of the error in the internal-layer subdomain 
were smaller in the computational results. Thus, the a priori error bounds of the 
asymptotic analysis did not reflect all of the accuracy and behavior of the computa- 
tional algorithm. 

Asymptotics identified two subdomains and provided preconditioners for the prob- 
lems within the subdomains. The preconditioner for the full equation in D IL is the use 
of the local scale x=(x-T)/e dictated by Theorem 2. This scale allows the diffusion 
to be modeled accurately, hence the grid is fine enough to resolve the shock. It is 
reasonable to use this scaling in the method, because computationally the internal- 
layer subdomain is of width O(e). The preconditioning in the outer-region subdomain 
D 0 r is to solve (8) in place of (2), and was justified in Theorem 2. First, the domain 
decomposition and preconditionings are combined with a functional iteration to form 
the computational method. The particular numerical methods discussed herein are 
not new; however, their combination to form this method is. 


5. DISCUSSION OF THE METHOD. An iteration is formed by linearizing 
the reduced problem. Each step of the iteration requires the solution of (8) in the 
outer-region subdomain, and (2) in the neighborhood of a shock. This is discussed in 
Section 5.1. Once the iteration has been described, the boundary detection scheme 
is presented. The method assumes no a priori information about the location of the 
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internal-layer boundary, and is supported by theory. In Section 5.3, convergence of 
the method is presented. Numerical details of the method will be presented in Section 

6 . 

5.1. Iteration. In general, each step of the iteration requires the solution of a 
linear convection-reaction equation in the outer-region subdomain followed by the 
solution of a nonlinear convection-diffusion-reaction equation in the internal layer. 
The convection-reaction equation 

( 17) u t fc+1 + U k U k+i - rU k+1 = o 

is formed by lagging the convection coefficient of (8). The boundary of the internal- 
layer subdomain is allowed to change during the iterations. Thus, denote the outer- 
region subdomain for iterate U k+l by D k OR , and denote the complement of D OR wi 
respect to D by D k IL . That is, U k+1 is obtained by solving (17) in D OR , then solving 

( 20 ) Thefsolution of (17) in D k OR is obtained via a modification of the method of char- 
acteristics to account for the forcing term rU k +K The characteristic transformation 
(x,t) — . (£, r) is defined by setting t = t, and by solving 


(18) 

with initial conditions 


if! = 

or 


x fc (0) = £, for b > i > 0; 
x fc (y a - a (0) = 0> for£<0; 

= 0, for £ > b; 

where (£ r) = (y a (r),r) is the image of the curve (x,t) = (0,i), and (£,r) — (yk( r )> T ) 
is the'image of the curve (x,t) = (M)- Utilizing this transformation, it is a simple 

task to solve 

(19) = rU k+1 

in place of (17) along the characteristics defined by (18). This transformation becomes 
singular in a neighborhood of a shock, hence cannot be applied in the internal-layer 
subdomain. This fact is the basis of the procedure used to determine dD JL . 

Solutions to the reduced equation are poor approximations to the solution of the 
full equation in regions of large gradients, such as in the internal-layer subdomain. 
Thus, the full equation is solved in the internal layer at each iteration. The equation 

(20) ’ U k+l + U k+1 U k+1 - eU k : k - rU k+1 = 0 

is solved in the internal-layer subdomain for each k. Boundary data for the eternal- 
layer subdomain is provided by the solution of (17) in the outer region. This is justified 
by observing that for c small enough, the boundary of the internal-layer subdomain 

will be an inflow boundary. 
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5.2. Boundary Detection. It is desirable to be able to compute the location 
of the internal-layer subdomain during the course of the iteration. This has the 
advantage of requiring less a prior information. In addition, if the initial guess provides 
a poor approximation to the location of the internal-layer subdomain, the method will 
be able to correct the location of the boundary in the course of the iteration. 

The method used to locate the internal-layer subdomain boundary is based on 
properties of the transformation used to solve (17). The transformation used to solve 
(17) will become singular (or nearly singular) in the region where characteristics merge 
(see Figure 2). Thus, the Jacobian 


(21) J k = dx/d $ , 

of the transformation (18) will be asymptotically small in a neighborhood of the shock, 
while it is O(l) in the outer-region subdomain. (For more details on the relationship 
between the magnitude of J k and the nature of the transformation, see [31].) This is 
the measure used to locate the boundary of the internal-layer subdomain. The size of 
the Jacobian is monitored via the solution of an ODE along each characteristic path 
of interest. Combining the partial of (18) with respect to £ and the partial of (21) 
with respect to r, the equation 


( 22 ) 


dj l 

dr 


= f 


dU k 

dx 




may be derived. The Jacobian is determined by solving this equation subject to 
•J*(*o><o) = J k for (x 0 , to) 6 II. The behavior of the solution inside the domain D is 
of primary interest; therefore, it is sufficient to monitor J k J k (x t ,ti)/j£ in place of 
J k to determine the boundary. The ratio is determined using (22) with J k in place of 
J , subject to J k (x 0 ,t 0 ) = 1 for (x 0 , to) € IT. The curve on which J k becomes nearly 
singular, that is, where 


( 23 ) J k {x,t) = TOL, 

for some small number TOL, is the boundary of the internal layer subdomain for 
iteration k. The subdomains separated by (23) will, in general, be different than 
the subdomains used in the theorems. Thus, the subdomains used in the numerical 
method are 


( 24 ) D k OR = {(x,t)|(x,f) e D, J k (x,t) > TOL}, 

and, 

( 25 > D) L = {(x,f)|(z,f) € D, j*(x,() < TOL). 

The theory applies provided 

< 26 ) d, l Cfl‘ l cj) 1 c...c £>,V c D k lL -, 
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however, convergence was observed when this relation failed, thus the constraint (26) 
is not a necessary condition for the computational method to converge. 

Heuristics, based on both accuracy and efficiency, are used to choose TOL. If 
TO L is too small, then accuracy will suffer. This is because the internal-layer subdo- 
main will be too small, and the data provided at the boundary of the internal-layer 
subdomain will have large perturbations as compared to the desired solution. If TOL 
is too large, the internal-layer subdomain will be too large, and the computational 
mesh will be refined in regions where the solution is smooth, creating excess work. 

5.3. Convergence. An advantage to this method is the availability of extensive 
error information. A global error bound based on Theorem 3 will be presented in this 
section. First, convergence of the iteration is established by showing the iteration is 
a contraction mapping. For more details on these results, see [32]. 

The convergence of the iteration (17) to a solution of (8) in the outer region will 
be shown by comparing successive iterates, then establishing a lower bound on the 
latest time at which the iteration is a contraction. For the sake of the theorem, the 
boundary of the internal-layer subdomain is assumed to be stationary from iteration 
to iteration {£>%£ = D k OR = D or ). With the analysis that follows, this theorem 
provides a lower bound for the largest time at which the iteration will converge. 

THEOREM 4. Let be the set of successive iterates of (17) in the sub- 

domain Dor satisfying (4-6) with initial guess U°. Assume U° satisfies (4-6) and is 
Lipschitz continuous on D. The boundary data are assumed to satisfy the compatiblity 
conditions (7), (9-10) and to have bounded first and second derivatives. Let 

6 = sup |Z7* — U k ~ l \. 

D 


Then 

(27) \U k+1 - U k \ < SCe- xt {e m - l) 

for ( x , t) e Dor . Here C, A and R are known positive constants. 

This theorem provides an upper bound on the latest time for which the iteration 
converges. Apply the infinity norm to (27) to obtain 

Then the following corollary provides the conditions for convergence. 

COROLLARY 5. Suppose that the conditions of Theorem 4 obtain. Let T max be 
the largest positive number such that 

C = sup Ce~ xt [e Rt — 1) < 1. 

<Tmfcx 
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If the bound on time in (S) satisfies 0 <T< T max , then the iteration in D or defined 
by (17) is a contraction mapping; therefore, the sequence of iterates converges to 
v = lim^oo U k = U°°, which is a solution of (8) on D 0 r satisfying (4-6). 

A statement of an a priori error bound for the computational method is presented 
in Corollary 6 below. As with Theorem 3, the bound is sharp in D 0 r\ however, the 
bound is crude in the region of the shock. 

COROLLARY 6. Let u be the solution to (2) satisfying (4-6). Suppose each 
iterate U k is obtained by first solving (17) in D or subject to (4-6), then solving (20) 
on D il with boundary data U k on dD IL . Assume that the compatibility conditions 
(7) > (9-1 0) obtain, and that the data (4-6) are bounded with bounded derivates. Suppose 
Q < T T maX j and let v — U°°. If E — ||u — v||i, then for e small enough 

E = O(e) 

xti Dor, and 


E = 0(e 1/4 ln 1/2 e) 


in Djl. Here A = 0(e*/ 4 ln 1 ^ 2 e). 

As with Theorem 3, the computational results reflect that both A and the error 
in the internal-layer subdomain are smaller than the corollary suggests. 


6. NUMERICAL DETAILS. The asymptotic analysis has provided a means 
to precondition the numerical problems. Because the sub-problems are well condi- 
tioned, the choice of numerical schemes may be made from a variety of standard 
methods. This is not usually the case. The class of problems for which the new 
algorithm is applicable are notoriously difficult to solve, and only a small number 
of schemes could be employed for its solution (prior to preconditioning). Since the 
sub-problems are well conditioned, numerical schemes used in the method presented 
here can be chosen based on criteria such as efficiency or the potential to exploit 
parallelism. 

6.1. Schemes in the Outer Region. The method of characteristics is used to 
solve the hyperbolic PDE in the outer-region subdomain. This method allows the 
exploitation of physically motivated parallelism. In addition, the method of charac- 
teristics allows handling of the free boundary at dD k IL in a straightforward manner. 

The method is not limited to using the method of characteristics for the problem 
* n D k L - For example, the schemes for hyperbolic conservation laws [1,8,35] might 
be modified to account for the term rU k and used. If this where done, then the 
method for detection of dD k IL could be based on the gradients instead of monitoring 
the Jacobian. 
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The method of characteristics scheme involves laying down a characteristic co- 
ordinate system, updating solution values, and monitoring the Jacobian. To update 
the solution and monitor the Jacobian requires a negligible computational cost. The 
discrete version of (18) is used to determine the characteristic coordinate system. 
Thus, for each iteration k, a set of characteristics {if },£i are computed and used 
as the computational grid. Here, the superscript k is the iteration number, and the 
subscript i identifies the characteristic. To determine a characteristic, 

(28) ^£ = P*(*?,<) 

is solved for each t = 1 to Ip. All of n is an inflow boundary for D, thus initial 
conditions 


if(r.) = £», 

are specified along all of n. On the (x,t) = (0, i) portion of II, the initial condition is 
= 0 at Ti = il, for * = 1 to I a . On the ( x,t ) = (z,0) portion, the initial condition 
is £, = (t - I a - l)h at Ti = 0, for i = I a + 1 to / 7 . And on the (z, t ) = (6, t) portion, 
the initial condition is & — b at r,- = (* — I~, — 1)/, for t = /-, + 1 to Ip. Here l and h 
are the increments for- time and space, respectively. The locations of characteristics 
are desired at time increments of At. It is assumed that r,- is some integer multiple of 
At] thus, the same temporal points are used for all characteristics for all iterations. 
Each characteristic is obtained on a At interval using the Trapezoid rule to solve (28). 
The Trapezoid rule is solved via a Newton Iteration. The computed value of zf(t n ) 
is denoted zf,-. 

As mentioned in Section 5.2, J k is monitored along each characteristic to deter- 
mine the boundary of the internal-layer subdomain at each iteration. The solution to 
(22) along characteristic t is 

(29) J fc (zf(t),t) = exp[5,^(t)], 
where 

(30) S,‘(i) = /'[7*(x‘(r),r)<(r. 

Since computations of (30) are better conditioned than those of (29), S k is monitored 
in place of J*. Denote the computed value of S k {tj) by S k j. The integral is determined 
using the right-hand rectangle rule 

(31) ^i,i + 1 = &i,j + Atf7*(zf ;+1 , iy+i) 

with the initial value S*(r,) = 0. The monitoring is performed at a minimal cost. 
It is necessary to keep only the most recent value of S k , thus minimal storage is 
required for this technique. In addition, the values of U k (x k j +1 , ty+i) are saved from 
the Newton iteration, hence very little computational cost is required. 
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The criteria used to determine the boundary will now be made more precise. The 
boundary dD) L is defined by J k = TOL, hence a characteristic is considered to be in 
the outer region as long as 

s,‘,- > ln(rOL). 

Let t = tj = jAt be the first time this inequality is violated for characteristic i 
during iteration k. Then the point (x, t) = (x*-.,j"Ai) is considered to be inside 

Dj L , and characteristic * is considered to be incident with the boundary at the point 
(x, t,) = (x^j. j, (j — l)At), and is flagged as being part of Dj L . 

After x k j has been determined, the solution U k+1 is computed by solving (19) 
subject to (4-6) using the right-hand rectangle rule 

(32) 0&*. « [1 + A( 

This formula is used until either j = T / At, or until characteristic i enters D) L , 
whichever happens first. This formula has minimal computational requirements; how- 
ever, the iteration requires the storage of the most recent iterate for a portion of D. 

6.2. Schemes In the Internal Layer Region. The sub-problem in the internal- 
layer subdomain requires the solution of a parabolic PDE subject to boundary data 
provided by the solution in the outer region. There are two major aspects of the 
computations in the internal-layer subdomain — mesh generation, and the difference 
technique. The mesh follows the shock, and has been scaled; therefore the variation in 
the solution is resolved on the new coordinate system. Thus, the mesh provides a pre- 
conditioning for the problem in the internal-layer, and the computations are not overly 
sensitive to the particular difference scheme used to solve the partial differential equa- 
tion. Russell’s Modified Method of Characteristics (MMC) [30], an explicit/implicit 
finite difference method, was chosen to solve the equation due to the regularity of the 
linear algebra problems which it generates. As with the schemes in the outer-region 
subdomain, other methods (see [7,5]) could be employed for the solution in b) L . It 
is necessary that the boundary of the internal-layer subdomain be identified with re- 
spect to the grid in the internal-layer subdomain. This is described first. Then the 
finite difference technique is reviewed. 

The base of the internal-layer subdomain is identified by finding which character- 
istic (or set of characteristics) has the lowest time at which it is flagged as being part 
of the internal-layer subdomain boundary. Denote the computed value of t r by t r . For 
t > t r , the base is taken to be the region between the two outer-most characteristics. 
This could result in non-flagged characteristics being part of the base, but this is not 
a problem. Once the base characteristics have been located, it is simple to identify 
whether a flagged characteristic is on the left or on the right boundary of DHj. A 
flagged characteristics with an index of lower value than a base characteristic is on 
the left boundary, and a flagged characteristics with an index of higher value than a 
base characteristic is on the right boundary. 
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A description of the method used to locate the left and right boundaries is fa- 
cilitated by first introducing the coordinate system. The computations will be one 
on a scaled and translated coordinate system with temporal variable t - t/e. (1 
spatial variable will be described later). The temporal grid is the set of points {<„}. 
where C = nAf for n = t r /(eA f) to T/(e Af). These points are a refinenwnto 
{!,}, the temporal points on which the characteristics are known. The left and right 
boundaries of the internal-layer subdomain at time t; are denoted by and 
respectively. Algorithm 1 describes the method used to determine fl; the procedure 
used to determine L is symmetric, and hence will not be described. 

Do n = i r /{eAV) to T/{e&t*) 

If (t; ^ tj/e) for any j 
then R n Rn-i 
otherwise 

if := number of characteristics incident with the right boundary at t n 

If H = 0 

then R n := R n - 1 
If H > 1 
then 

i := index of right-most characteristic incident 
with the right boundary 

R n == x i,j 

ALGORITHM 1 Determination of R n . 

It is now appropriate to describe the coordinate system on which the computations 
in the internal layer are based. Denote the middle of the internal-layer subdoxmin a 
time t by M(t) = [R{t) + L{t))/ 2. Then the spatial coordinate in the internal layer 

(33) x* = [x - M(t)]/e. 

The temporal coordinate is also scaled, t* = t/e. Equation (20), may be written as 


(34) 


Uf. +1 + 


fdM 
V dt 


u k+1 ) U k J x - u£\ - efU k+i = o 


Here, U k+ 1 {x\f) = U{M + ex\et *), and r(i*) = r(M + ex*). This is the form of the 

equation solved in the internal layer. . . tVip rpader 

Equation (34) on the grid defined above is solved using the MMC. The reader 

should refer to [27,29,30] for a more complete description of the MMC; how ' ver ' “ 
review of the method is presented here for the sake of completeness. Denote tlmthe 
spatial points of the computational grid by x, - »Ax , for i .... 

Ax‘ is a constant and the width of D) l varies with time, the number of grid p 
also varies with time, and /“ = /“(•). The MMC involves approximation of the 
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convective term 0*» + (* - &‘«) by a backward Euler approximation along 
the subcharacteristics of (34) using the following formula 


(35) 


£tf +1 (xV*) + - tf t+1 (xV)j t£ +1 (x*,f*) 

- (^* +1 (x*,r) - U k+ 1 (z*,t* - A**)]/At*. 


Here + At*(^£^n - - At*)) and **«(*•,!• - At*) are 

e ermined by linear interpolation between spatial grid points. This linear interpo- 
lation is performed element by element, thus there is ample opportunity to exploit 
parallelism within the method here. Once this quantity has been calculated, the s< 

ond derivative is approximated using a centered difference. Hence, the full formula 
the internal layer is 


sec- 
in 


(36) 


1 - cA t*f{x*) + 


(Ax*) 


-At* 


U k+1 (x*,t*)~ 


At 


(Ax*) 2 ( Uk+1 ( x * + Ax*,r) + U k+1 (x* - Ax* , t*)) = U k+1 (i*,t* - At*). 

There is a domain of dependency requirement on the method [29], which requires the 
absolute value of the partial with respect to x* of the convection coefficient to be 
bounded by 1/At*. Since M is independent of x*, this requirement is that 

< 37 ) \U£ l {x*,t*)\<l. 

Extra inner iterations may be necessary to step the solution between the temporal 
grid lines, £*. * 


7 . OUTLINE AND IMPLEMENTATION OF THE ALGORITHM. 

o u J’ 1 ' Algonthm - As a summary, the numerical method is outlined in Algorithm 
e ow. e parameters such as TOL and At are assumed to have been provided by 
the user. Several steps are not mentioned. For example, if the domain of dependency 
requirement is violated in the MMC, it may be necessary to reset At*. However 
Algorithm 2 shows the major computational requirements. 

The algorithm requires an initial guess. In the theory, the initial guess must 
satisfy the boundary conditions (4-6), and must be Lipschitz continuous; however 
for the computation^ method, it is only required that some approximation technique 
which is consistent with the effects of the terms of (2) be used to determine U° The 
MMC was chosen for the initial guess. The computations were done on the original 
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Do k = 1 till converged 

I. Solve in D k OR 

A. Do i = 1 to Ip 

1. j := Ti/At 

2. x k ■ £,• 

3. Do while j < T / At and S k < ln(TOL) 

a. j := j + 1 

b. Step x k from tj- \ to tj 

c. Compute solution value using (32) 

d. Update size of S k with (31) 

II. Solve in b) L 

A. Determine dD) h 

1. Find base (determine t r ) 

2. Find R using Algorithm 1 

3. Find L (symmetric to iZ) 

B. Set boundary values along dD k L 

1. Initial conditions 

2. Right boundary values 

3. Left boundary values 

C. Discretize 

1. t* := t r /e 

2. Determine At* to satisfy (37) as a subdivision of {tj} 

3. t* := t* + At* 

4. Do while t* < T/e 

a. t* := t* + At* 

b. Obtain U k at t* 

i. Form right hand side of (36) 

ii. Solve implicit portion of (36) 

c. If (37) is violated, go to step C.l. 

ALGORITHM 2 Computational method. 
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coordinates, x, and t. In order that the domain of dependency requirement be met, 
artificial diffusion was added instead of restricting the size of At. 

Use of artificial diffusion could lead to some ill effects, especially in the location 
of the shock. The size of the diffusion coefficient effects the speed of a nonsteady 
shock. Thus, using artificial diffusion may result in computing the location of a 
nonsteady shock incorrectly. This will be shown in the Model Problem II. However, 
this inaccuracy is acceptable because the method is not sensitive to errors in the initial 
guess, as long as U° is continuous. 

7.2. Parallelism. Parallelism may be exploited at several levels in the imple- 
mentation of Algorithm 2. Consider first the parallelism which may be exploited 
in the solves for the characteristics. Characteristic xf is obtained at discrete points 
in time by solving (18) using a Newton iteration. These solves may be scheduled 
asynchronously, or they may be grouped as vectors. Grouping characteristics and 
assigning the spatial location of each characteristic in a group to a component of a 
vector allows the exploitation of vector processing capabilities. Once the location of 
the characteristic is known for a time step, the value of ln(J*) may be approximated 
using (31). 

For a large number of processors, a parent process could spawn a task for each 
characteristic, and allow them to all execute in parallel. In turn, the task solving a 
particular characteristic could then spawn two tasks, one for the monitoring of the 
Jacobian, and one for the updating of U k+1 . To avoid extra computations, the child 
task monitoring the Jacobian would need a means of interrupting the parent task 
computing the characteristic; however, no communication from the task computing 
U k+1 to the parent task is needed. This is reflected by the data dependency graph in 
Figure 3. The parent routine is labeled OR. The tasks which compute characteristics 
are labeled C„ for t = 1 to I p . Each characteristic task spawns two processes, one for 
the Jacobian monitor, and one to obtain the value of the solution. These are labeled 
J and U, respectively. The data dependency graph for D^ R is a subgraph of the 
dependency graph for the whole problem. 

For a smaller number of processors, since a large number of characteristics will 
be computed, the exploitation of the medium grain parallelism outlined above is not 
needed. Thus, for the implementation on the Alliant FX/8, a single task performs the 
Newton iteration to determine the new location of xf. Then the same task updates 
the values of S k and U k+l at the new characteristic location. 

Other parallelism evident from the description of Algorithm 2 is reflected in Figure 
4. Nodes are labeled with the corresponding step number of Algorithm 2. The only 
sequential step is the location of the base of dD k IL , and this is a small portion of the 
overall computations. The major portion of the computations are the characteristic 
solves in the outer region and the discretization for the internal layer. More efficient 
methods could possibly be applied in the internal layer, potentially providing more 
parallelism. 

A less obvious source of parallelism is the exploitation of another type of domain 
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FIG. 3. Data dependency for characteristic solves. (From [32].) 

decomposition to form a pipeline out of the outer (fc) iteration. To obtain U ! up 
to time t = T, requires the knowledge of U' for 0 < i < T„ but not for t > T,. 
Thus, U 1 for Ti < t < T? can be computed at the same time that U is determine 
for 0 < t < Tj. In general, subdivide the domain D in the temporal direction as 
0 < Ti < T 2 . . . < Tff < T. While U 1 is being computed for T k < t < T k + 1, the 
computations for U 2 through U k could be taking place, thus forming a pipeline based 
on the temporal subdivision of the domain. 

7.3. Implementation. Algorithm 2 was implemented on an Alliant FX/8 using 
a package called Schedule [14]. This package provides a common user interface to the 
parallel capabilities of a variety of shared memory parallel computers All of the 
synchronization required to enforce the data dependencies is automatically provided 
by the Schedule Package once the graph has been specified correctly. Moreover, 
there are no machine dependent statements within the user code. All such machine 
dependencies are internal to Schedule. This provides for transportability of the code 
between the various machines Schedule has been ported to. The implementation is 
meant as a demonstration of the viability of the method. Thus, not all of the available 
parallelism has been exploited. Even so, significant speedups were achieved, and will 

be discussed in the section on the experiments. 

A useful feature of Schedule is the automatic generator of a data flow graph 
associated with a computation. In Figure 5, the graph Schedule produced for an 
older version of the code [23] is shown. The nodes represented by circles form the 
static portion of the data dependency graph (the portion of the graph which does no 
change for different input data). The rectangular nodes associated with static nodes 
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FlG. 4, Data dependency for algorithm. 
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FlG. 5. Dependency graph generated by Schedule . (From [23].) 


represent dynamically spawned processes which, in this case, are each a characteristic 
grid line. The graph represents a snapshot of computation shown midway through 
its execution. The black nodes represent computed processes, the hatched nodes 
represent active processes, and the white nodes represent processes waiting to execute. 
Information from this graph and various statistics available with it have been used to 
improve load balance in this algorithm. 

The matrix equation in step II.C.4.b.ii is a symmetric positive definite tridiagonal 
matrix of about 300 unknowns. Since the number of unknowns is small, the Linpack 
tridiagonal solver, sptsl, [15] was parallelized instead of using a more complicated 
block scheme. The scheme in sptsl uses two dual sweeps of the tridiagonal matrix. 
The first dual sweep is for the forward elimination, and begins at both the top and 
bottom of the matrix, then works inward. The backward substitution sweep begins 
at the middle of the matrix, then works outward in both directions. Thus, the com- 
putations can be parallelized for a maximum speedup of two. The implementation 
used for this method had a speedup of about 1.5 on 2 to 8 Computational Elements 
(processors) of the FX/8 over the compiler-optimized version of the sequential code. 


8. EXPERIMENTS. Two problems are solved, each demonstrating different 
features of the algorithm. Model Problem I has a steady shock and no forcing term, 
with an exact solution being known. It demonstrates that the behavior of the error in 
the computations is the same as the theoretical error estimates of the theorems in the 
outer regions. Model Problem II has a nonsteady shock. It is used to demonstrate 


that the method is not sensitive to the initial guess. 

8.1. Model Problem I. The first model problem will demonstrate that the 
error tends to zero as e j 0. With no forcing term, (2) is Burgers’ equation, 

(38) t if t uUj — 0* 

The initial and boundary conditions are that 

u = —2 tanh(x/c), 

for (x, t) £ IT. Under these conditions, the exact solution to (38) is u = — 2tanh(x/c), 
hence the computed solution may be compared easily to the exact. The Li error is 
presented for runs with different values of e in Table 1. The error was measured at 

Table l 2 Li error 


6 

Error 

h- 1 

o 

1 

to 

5.0 x 10" 02 

10" 3 

2.2 x 10-° 3 

10" 4 

2.3 x 10~ 04 


time t = .1, and remained constant for the remainder of the computational domain. 
Only one pass of the domain decomposition was needed for this problem. The results 
presented in Table 1 indicate that the error in the computational method is O(e). 

8.2. Model Problem II. The second model problem has a nonsteady shock, 
and shows the effect of the forcing term on the location of the shock. In addition, 
this method demonstrates that the method is not overly sensitive to the accuracy of 
the initial guess. The equation solved is 


u t + uu x — eu xx — 4 sin(2?rx)u = 0, 


subject to the initial guess 

u(x, 0) = cos(ttx), for 0 < x < 1 = b 
(see Figure 6). Then the boundary conditions are 

u(0, t) = 1, u(l,t) = -1, for T > t > 0. 

For the experiments presented here, c = .001. If there were no forcing term (r = 0), 
then this problem would have a steady shock develop in the center of the domain; 
however, the forcing' term is r = 4 sin(27rx) . This represents a duct of width A(x) = 
exp[2 cos(27rx)/7r], which has the general shape of the one in Figure 1. The effects of 
the shape of the duct on the location of the shock are reflected by the internal-layer 
subdomain having a different location as time increases. The movement of the internal 
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layer subdomain can be seen in Figure 7, where the characteristics of the outer-region 

solution and internal-layer subdomain after four passes of the domain decomposition 
are shown. 

The initial iterate is obtained using the Modified Method of Characteristics on 
a rectangular grid with 100 spatial points. Figure 8 shows the initial iterate, U°, at 



* 1 * _L_ 

0.0 0.5 1.0 


X 

Fig. 8. Initial guess at t = .3. (From [32].) 

time t = .3. The use of artificial diffusion resulted in the location and the width of 
the shock being wrong for the initial iterate. Comparing the initial iterate, U°, with 
U in Figure 9, these errors may be seen. The initial guess has the wrong amplitude, 
and the shock is slightly to the left and wider than the shock in U 3 . 

The succession of internal-layer subdomains may be seen in Figures 10-13. 

One of the manifestations of the convergence of the iteration is that the internal- 
layer subdomain has a much smoother boundary after convergence than before. The 
boundary for U* is smooth up to approximatedly time t = . 3 . 

Timings for computing the initial guess and the first iteration are presented in 
Table 2. Speedup is the execution time for multiple CEs divided by the time for one 
CE. Efficiency is the speedup divided by the number of CEs. The data indicates that 
significant speedups were attained, even though a significant portion of the parallelism 
was not exploited. For example, the pipelining described in Section 7.2 was not used. 
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9. CONCLUSION. In this paper, asymptotics and numerics have been blended 
to form a new computational method. The method has potential to exploit a large 
amount of parallelism and provides high accuracy. Asymptotic analysis provided a 
theoretical basis for the domain decomposition, and guided in the derivation of rigor- 
ous local and global error bounds. 

Two types of subdomains were identified by the asymptotic analysis of this prob- 
em: smooth outer regions, and an internal-layer subdomain with a shock. Pipelining 
the outer iteration provided large-grain parallelism. Smaller-grain parallelism was ex- 
ploited by using a modification of the method of characteristics in the outer regions 
and by blocking the computations in the shock layer. 

The method was developed for an important problem in computational fluid dy- 
namics; however, the method is suitable for a wide range of problems in physics and 
chemistry. Namely, this approach is suitable for problems with internal layers and 
boundary layers interspersed with regions where the solution is smooth. Examples of 
such problems other than transonic flow through a duct include the location of stagna- 
tion points (flow of zero velocity) where two opposing jets intersect, and combustion 
ironts. 

The availability of estimates and bounds on the error is important in the design of 
numerical methods. Rigorous a priori error bounds were established for the method 
presented here. In addition, the particular numerical schemes used for the subprob- 
lems allowed a posteriori error estimation. The a priori error bounds were shown to 
much larger than the errors observed in the computations; thus, sharper error bounds 
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NOTATION 


D 

(3) 

Dil 

(14) 

Dor 

(15) 

Dil 

(25) 

Dor 

(24) 

L(t) 

Sec. 6.2 

Ln 

Sec. 6.2 

p 

(2) 

p b 

(8) 

m 

Sec. 6.2 

Rn 

Sec. 6.2 

sm 

(30) 

sti 

(31) 

T 

(3) 

U k 

(17), (20) 

t r 

Sec. 4 

t* 

Sec. 6.2 

u 

(2) 

u 

(8) 

X 

Sec. 4 

X* 

(33) 

r 

(35) 

x k (t) 

(18) 


Sec (6.1) 

x t,J 

Sec (6.1) 


Computational domain. 

Internal-layer subdomain. 

Outer-region subdomain. 

Computational internal-layer subdomain. 

Computational outer-region subdomain. 

Location of the left boundary of D r L. 

Computed location of L(t) at t n . 

The full nonlinear operator. 

The reduced nonlinear operator. 

Location of the right boundary of DjL. 

Computed location of R{t) at t n . 

Integral for monitoring Jacobian. 

Discrete values of S k {t). 

Upper bound on time t for D. 

Time that the solution to the reduced problem becomes multivalued. 
Scaled and translated internal-layer coordinate. 

The solution to the full operator. 

The solution to the reduced operator. 

Scaled and translated internal-layer coordinate. 

Spatial coordinate for internal-layer computations. 

Spatial location used in MMC. 

Characteristic coordinate. 

Characteristic grid line i for iteration k. 

Computed value of x k at time in- 
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